Selection among males reduces mutation load

Load packages

library(tidyverse) # tidy style coding
library(brms) # Bayesian models
library(bayesplot) # pretty bayes visuals
library(tidybayes) # Bayesian aesthetics
library(loo) # to use information criteria in brms models
library(MetBrewer) # colours
library(rcartocolor) # more colours
library(pander) # tables
library(kableExtra) # tables
library(patchwork) # putting plots together
library(DT) # for search- and saveable tables
library(ggdist) # for ribbon plot
library(png) # to load images
library(grid) # to plot images

Supplementary methods

\(~\)

Table S1. Recipe for food medium used in our experiment. The provided quantities make ~ 1 litre of food.

tibble("Ingredients" = c("Soy flour", "Cornmeal", "Yeast", "Dextrose", "Agar", "Water", "Tegosept", "Acid mix (4 mL orthophosphoric acid, 41 mL propionic acid, 55 mL water to make 100 mL)"),
       "Quantity" = c("20 g", "73 g", "35 g", "75 g", "6 g", "1000 mL", "17 mL", "14 mL")) %>% 
  pander(split.cell = 40, split.table = Inf)

Ingredients	Quantity
Soy flour	20 g
Cornmeal	73 g
Yeast	35 g
Dextrose	75 g
Agar	6 g
Water	1000 mL
Tegosept	17 mL
Acid mix (4 mL orthophosphoric acid, 41 mL propionic acid, 55 mL water to make 100 mL)	14 mL

\(~\)

img <- readPNG("Figure_S1.png")
 
grid.raster(img)

Figure S1. Crossing scheme used to integrate the GFP constructs and ap^XA marked translocated second and third chromosome balancers into the LH_M genetic background. We replicated the crosses 12 times to supply the flies used in generation zero of experimental evolution; 6 times using the Ubi GFP construct and 6 times with the 3xP construct. G = generation.

\(~\)

Analysis

\(~\)

Mortality analysis

\(~\)

We conducted single-pair full-sibling crosses to erode heterozygosity across the genome, thereby exposing mutation load. Our proxies for mutation load - extinction rate and productivity decline - may also have been affected by differences in the intensity of interlocus sexual conflict between subpopulations carrying autosomes with different selection response histories. We specifically designed our experiment to minimise this effect, by enforcing monogamy between all breeding individuals across the extinction assay. However, we cannot discount males carrying male-limited autosomes and/ or control autosomes being more harmful to females than males carrying autosomes with a female-limited selection response history. This might occur because male harm is a by-product of male-male competition, which cannot respond to selection in the female-limited treatments. Furthermore, traits genetically correlated with harm such as activity levels or metabolic rate may be selected in the opposite direction among females, as these traits are highly sexually dimorphic. This may drive the evolution of reduced harm during female-limited evolution.

We observe 135 female (6.2%) and 21 (1%) male mortalities across the 2184 Drosophila vials used during the 20 generations of the experiment (we do not include the 9 final generations observed in Block 1 in this analysis). Male mortality did not occur often enough to provide the power required for formal analysis. We instead focus on modelling female mortality and test whether this is affected by selection response history.

We test this hypothesis to assess whether it is required to censor subpopulations where the breeding female died in the generation of extinction.

\(~\)

Load in the data

mortality_data <- 
  read_csv("Data/Mortality_data.csv") %>% 
  rowid_to_column("subpopulation") %>% 
  pivot_longer(cols = 8:27, names_to = "Generation", values_to = "Female_mortality") %>% 
  mutate(
    Generation = as.integer(str_remove(Generation, "Gen_")),
    Week = case_when(
    Block == 1 ~ Generation,
    Block == 2 ~ as.integer(Generation + 9)),
    Week = as.factor(Week),
    across(1:7, as.factor),
    Male_mortality = Female_mortality,
    Female_mortality = if_else(Female_mortality == "FEMALE" | Female_mortality == "BOTH", 1, 0),
    Male_mortality = if_else(Male_mortality == "MALE" | Male_mortality == "BOTH", 1, 0)) %>% 
  rename(Evolution_treatment = Treatment)


# Create a function to build HTML searchable tables

my_data_table <- function(df){
  datatable(
    df, rownames=FALSE,
    autoHideNavigation = TRUE,
    extensions = c("Scroller",  "Buttons"),
    options = list(
      dom = 'Bfrtip',
      deferRender=TRUE,
      scrollX=TRUE, scrollY=400,
      scrollCollapse=TRUE,
      buttons =
        list('pageLength', 'colvis', 'csv', list(
          extend = 'pdf',
          pageSize = 'A4',
          orientation = 'landscape',
          filename = 'mortality_data')),
      pageLength = 8500
    )
  )
}

my_data_table(mortality_data %>% 
                select(Mother_strain, Father_strain, Cross, subpopulation, Block, Week, Generation, Evolution_treatment, Female_mortality, Male_mortality))

\(~\)

Column explanations

Mother strain: we measured extinction and productivity of families derived from 36 different crosses. Crosses were only conducted between strains that experienced the same Evolution_treatment. This column shows the strain that the founding female of the subpopulation was derived from.

Father strain: we measured extinction and productivity of families derived from 36 different crosses. Crosses were only conducted between strains that experienced the same evolution treatment. This column shows the strain that the founding male of the subpopulation was derived from.

Cross: an identifying ID for each of the 36 combinations of Mother strain and Father strain.

subpopulation: lines of flies founded by single inseminated females in generation zero. We generated the next generation of each subpopulation by crossing a single full-sibling dyad.

Block: the experiment was run in 2 distinct blocks, using flies separated by 9 generations.

Week: the week that the generation of each block was conducted in. Laboratory conditions may have differed slightly between weeks. For example, a new batch of food was cooked for each week. Note that blocks overlap between weeks 9-20.

Generation: the generation of the extinction assay, ranging from 1 to 20.

Evolution_treatment: the strains had been exposed to one of three evolutionary conditions for 20 generations: a female-limited response to selection, a male-limited response and a control condition where an evolutionary response occurred in both sexes.

Female_mortality: a 0 indicates the female survived until she was removed from the 3-day breeding vial. A 1 indicates she died during this period. Females were 1-4 days old when they entered the breeding vial.

Male_mortality: see Female_mortality, except this column documents male mortality in the breeding vials.

\(~\)

Modelling approach

\(~\)

We fit a binomial model with Female_mortality as the response, Evolution_treatment as a fixed effect and Week, subpopulation and Cross as random effects.

Priors

We fit moderately informative priors, with the aim to regularise our posterior estimates and improve model fitting by ruling out non-sensical values.

Note that each of these priors is expressed on the logit scale, due to the binomial logit link.

\(\alpha\) (the intercepts) ~ Normal(\(\mu\) = -2, \(\sigma\) = 2)

\(\beta\) (the fixed effects) ~ Normal(\(\mu\) = 0, \(\sigma\) = 2)

\(\sigma\) ~ Exponential(rate = 1)

female_mortality_model <- 
  brm(Female_mortality ~ 1 + Evolution_treatment + (1|Week) + (1|Cross) + (1|subpopulation),
      data = mortality_data %>% filter(Female_mortality != "NA"),
      family = bernoulli,
      prior = c(prior(normal(-2, 2), class = Intercept),
                prior(normal(0, 2), class = b),
                prior(exponential(1), class = sd)),
      iter = 8000, warmup = 4000, chains = 4, cores = 4, seed = 1,
      control = list(adapt_delta = 0.98),
      file = "Fits/female_mortality_model")

Extract model predictions and plot

mortality_predictions <-
  female_mortality_model %>% 
  as_draws_df() %>% 
  mutate(Control = inv_logit_scaled(b_Intercept)*100,
         `Female-limited` = inv_logit_scaled(b_Intercept + b_Evolution_treatmentFemale)*100,
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Evolution_treatmentMale)*100) %>% 
  select(Control, `Female-limited`, `Male-limited`)
    
# put in easy to plot format

mortality_predictions_long <-
  mortality_predictions %>% 
  pivot_longer(cols = everything(), names_to = "Evolution_treatment", values_to = "Mortality")

# calculate the differences between each treatment

mortality_diff <-
  mortality_predictions %>% 
  mutate(`Control - Male-limited` = Control - `Male-limited`,
         `Male-limited - Female-limited` = `Male-limited` - `Female-limited`,
         `Control - Female-limited` = Control - `Female-limited`) %>% 
  pivot_longer(cols = 4:6, names_to = "Difference_contrast", values_to = "Difference") %>% 
  select(contains("Diff"))

Build Figure S2

mortality_plot <-
  mortality_predictions_long %>% 
  ggplot(aes(x = Mortality, y = Evolution_treatment)) +
    stat_halfeye(aes(fill = Evolution_treatment), .width = c(0.66, 0.95), alpha = 1,
               point_interval = "median_qi", point_fill = "white", 
               shape = 21, point_size = 4, stroke = 1.5) +
  scale_fill_manual(values = met.brewer("Hiroshige", 5)) +
  labs(x= "Female mortality (%)", y = "Inheritance treatment") +
  theme_bw() + 
  theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
        plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
        panel.grid.minor.x = element_blank(),
        legend.position = "none", #transparent legend panel
        text = element_text(size=14))

mortality_diff_plot <- 
  mortality_diff %>% 
  ggplot(aes(x = Difference, y = fct_relevel(Difference_contrast, "Control - Female-limited", "Male-limited - Female-limited", "Control - Male-limited"))) +
  stat_halfeye(aes(fill = Difference_contrast), .width = c(0.66, 0.95), alpha = 1,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.5) + 
  scale_fill_manual(values = c(carto_pal(7, "Peach")[2], carto_pal(7, "Purp")[1], carto_pal(7, "TealGrn")[1])) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", linewidth = 1) +
  labs(x = "Diff. in mortality (% points)", y = "Treatment contrast") +
  scale_x_continuous(breaks=seq(-2, 6, 2), limits = c(-3, 7)) +
  theme_bw() + 
  theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
        panel.grid.minor.x = element_blank(),
        plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
        legend.position = "none", #transparent legend panel
        text = element_text(size=14))

  mortality_plot / mortality_diff_plot +
  plot_annotation(tag_levels = 'a')

Figure S2. Female mortality is less frequent in subpopulations with female-limited selection response histories, suggesting that male harm may be less intense in these subpopulations. Panel a shows the posterior distribution of the mean percentage of female mortality events across the total number of vials that housed subpopulations throughout the extinction assay, split by selection response history. Panel b shows the posterior distribution of the difference between each treatment. The points show the estimated median, with associated 66 and 95% credible intervals.

\(~\)

To avoid this confounding our measure of mutation load, we can censor extinction events that co-occur with female mortality.

\(~\)

Extinction analysis

Load in the extinction data

\(~\)

extinction_data <- 
  read_csv("Data/Extinction_data.csv") 

# Block 1 runs for 29 generations, while Block 2 only runs for 20. To calculate the censoring variable, we need to split these by Block, mutate the data, then rebind them 

Block_1 <-
  extinction_data %>% 
  filter(Block == "1") %>% 
  # here we create a censoring column. If the family a) escaped or was killed by something unrelated to the experiment or b) survived the 20 generations of the experiment, then we code a value of 1. If the family went extinct, we code a value of 0. This allows us to right censor the data, thereby preserving the information it provides on extinction.
  mutate(across(Gen_1:Gen_29, ~replace_na(.x, "Escape")),
         Censored_alive = if_else(Gen_29 == "YES", 1, 0),
         Censored_escape = if_else(Gen_29 == "Escape", 1, 0),
         Censored = Censored_alive + Censored_escape,
         across(Gen_1:Gen_29, ~if_else(.x == "YES", 1, 0)))
        

Block_2 <-
  extinction_data %>% 
  filter(Block == "2") %>% 
  # here we create a censoring column. If the family a) escaped or was killed by something unrelated to the experiment or b) survived the 20 generations of the experiment, then we code a value of 1. If the family went extinct, we code a value of 0. This allows us to right censor the data, thereby preserving the information it provides on extinction.
  mutate(across(Gen_1:Gen_20, ~replace_na(.x, "Escape")),
         across(Gen_21:Gen_29, ~replace_na(.x, "Not measured")),
         Censored_alive = if_else(Gen_20 == "YES", 1, 0),
         Censored_escape = if_else(Gen_20 == "Escape", 1, 0),
         Censored = Censored_alive + Censored_escape,
         across(Gen_1:Gen_29, ~if_else(.x == "YES", 1, 0)))

# combine the Blocked data back into a single tibble

extinction_data_wrangled <-
  rbind(Block_1, Block_2) %>%
  mutate(across(1:7, as.factor), 
         Gens_to_extinct = Gen_1 + Gen_2 + Gen_3 + Gen_4 + 
           Gen_5 + Gen_6 + Gen_7 + Gen_8 + Gen_9 + Gen_10 + 
           Gen_11 + Gen_12 + Gen_13 + Gen_14 + Gen_15 + Gen_16 + 
           Gen_17 + Gen_18 + Gen_19 + Gen_20 + Gen_21 + Gen_22 + 
           Gen_23 + Gen_24 + Gen_25 + Gen_26 + Gen_27 + Gen_28 + Gen_29 + 1) %>% 
  rename(Evolution_treatment = Treatment, subpopulation  = ID) %>% 
  select(Mother_strain, Father_strain, Cross, subpopulation, Block, 
         Evolution_treatment, Gens_to_extinct, Gen_1:Gen_29, Censored_alive, 
         Censored_escape, Censored)

# Find the extinctions that co-occur with female mortality

extinction_mortality <-
  left_join(
    extinction_data_wrangled %>% 
      pivot_longer(cols = 8:36, names_to = "Generation", values_to = "Extant") %>% 
      mutate(Generation = as.integer(str_remove(Generation, "Gen_"))),
    
    mortality_data
  ) %>% 
  filter(Extant == 0 & Female_mortality == 1) %>% 
  mutate(Censored_mortality = 1) %>% 
  select(subpopulation, Censored_mortality)

extinction_data_wrangled <- left_join(extinction_data_wrangled, extinction_mortality) %>% 
  mutate(Censored_mortality = if_else(is.na(Censored_mortality), 0, 1),
         Censored_2 = Censored + Censored_mortality)

# Create a function to build HTML searchable tables

my_data_table <- function(df){
  datatable(
    df, rownames=FALSE,
    autoHideNavigation = TRUE,
    extensions = c("Scroller",  "Buttons"),
    options = list(
      dom = 'Bfrtip',
      deferRender=TRUE,
      scrollX=TRUE, scrollY=400,
      scrollCollapse=TRUE,
      buttons =
        list('pageLength', 'colvis', 'csv', list(
          extend = 'pdf',
          pageSize = 'A4',
          orientation = 'landscape',
          filename = 'extinction_data')),
      pageLength = 500
    )
  )
}

my_data_table(extinction_data_wrangled)

Column explanations

Mother strain: we measured extinction and productivity of families derived from 36 different crosses. Crosses were only conducted between strains that experienced the same Evolution_treatment. This column shows the strain that the founding female of the subpopulation was derived from.

Father strain: we measured extinction and productivity of families derived from 36 different crosses. Crosses were only conducted between strains that experienced the same evolution treatment. This column shows the strain that the founding male of the subpopulation was derived from.

Cross: an identifying ID for each of the 36 combinations of Mother strain and Father strain.

subpopulation: lines of flies founded by single inseminated females in generation zero. We generated the next generation of each subpopulation by crossing a single full-sibling dyad.

Block: the experiment was run in 2 distinct blocks, using flies separated by 9 generations. Block 1 ran for 29 generations, whereas Block 2 ran for 20 generations.

Evolution_treatment: the strains had been exposed to one of three evolutionary conditions for 20-29 generations: a female-limited response to selection, a male-limited response and a control condition where an evolutionary response occurred in both sexes.

Gens_to_extinction: the number of generations the subpopulation survived until it went extinct.

Gen_1:29: a 1 indicates the subpopulation was extant for the generation in question, while a 0 indicates extinction. Note that we continued propagating surviving subpopulations from Block 1 until generation 29, which coincided with generation 20 for subpopulations belonging to the second Block.

Censored: if the subpopulation 1) ended because flies escaped or were killed by something unrelated to the experiment, 2) went extinct because of breeding female mortality or 3) was extant after the final generation of the experiment, then we coded a value of 1. If the subpopulation went extinct, we coded a value of 0. This allows us to right censor the data as per the brms syntax, thereby preserving the information these censored cases provide on extinction.

\(~\)

Modelling approach

\(~\)

We fit a weibull survival model to estimate , the extinction rate of the subpopulations used in our experiment. The weibull model is an extension of an exponential decay model, and includes an additional shape parameter, which allows to vary across generations of inbreeding.

We right censor subpopulations that 1) ended because flies escaped or were killed by something unrelated to the experiment, 2) went extinct because of breeding female mortality or 3) were extant after the final generation of the experiment

Response variable

Gens_to_extinct is our response - a continuous variable that runs from 1 to 30 generations. Note that we measured extinction up until generation 29 in Block 1 and generation 20 in Block 2. A value of 30 indicates that the subpopulation was extant at the end of the experiment in Block 1, while a value of 21 indicates this in Block 2. We also include right censoring, to account for subpopulations extant at the end of the experiment, or that were removed from the experiment by a handling error. Censoring enables the inclusion of this ‘incomplete’ data in the model.

Fixed effects

Evolution_treatment: we include this to test for a causal effect of sex-limited experimental evolution on mutation load.

Block: extinction rate might differ between the blocks we split our experiment up into e.g. because of microvariation in the lab brought about because of a change of season, slight inconsistencies in Drosophila food consistency etc.

Varying/Random effects

Cross: we initiated each subpopulation by crossing two strains belonging to the same evolution treatment. There were 12 levels of cross nested within evolution treatment, for a total of 36 crosses. We kept subpopulations from the same Cross and Block in the same column of the Drosophila vial trays, so incuding this variable also controls for micro-environmental variation caused by within-tray position. We used this random effect in our first model.

Mother_strain and Father_strain: Used as an alternative to Cross in our second model. In this case, we fit a multi-membership model, where each subpopulation simultaneously belongs to two levels of the Cross random effect. That is, each subpopulation was founded my a mother from \(strain_i\) and a father from \(strain_j\) and they therefore always belong to two levels of strain.

Multi-membership vs single-membership

While the multi-membership approach may initially seem the more powerful of the two random effect structures, we favour the single-membership approach for several reasons. First, our primary aim is to estimate the overall causal effect of evolution treatment on generations to extinction, rather than that of each strain (nested within evolution treatment). Given that we balanced our crossing design so that each strain is equally represented in the experiment, estimation for each strain is not integral. Second, Cross captures nuisance environmental variation that Strain does not. We distributed subpopulations throughout three Drosophila vial trays by placing subpopulations from the same cross in the same column. Across the columns we split subpopulations up by evolution treatment, such that column one contained subpopulations from the female-limited treatment, column two contained subpopulations from the male-limited treatment and column three contained subpopulations from the control treatment. We repeated this pattern until all 36 crosses filled a column. As a strain is used in multiple crosses, it is widely distirbuted throughout the trays and therefore does not capture a location effect like cross does. We are of the opinion that controlling for this environmental variation is of greater importance than strain level differences in mutation load, which both cross and evolution treatment also contain information for.

Finally, we fit both models and use leave one out (LOO) cross validation to formally assess which model has better expected predictive accuracy, both in and out of sample.

Priors

We fit moderately informative priors, with the aim to regularise our posterior estimates and improve model fitting by ruling out non-sensical values.

Note that each of these priors is expressed on the log scale, due to the weibull log link.

\(\alpha\) (the intercepts) ~ Gamma(shape = 5, rate = 3)

This intercept prior predicts an appropriate starting point for per generation extinction rate on the log scale.

\(\beta\) (the fixed effects) ~ Normal(\(\mu\) = 0, \(\sigma\) = 1)

\(\sigma\) ~ Exponential(rate = 1)

\(k\) (The weibull shape parameter) ~ Normal(\(\mu\) = 0, \(\sigma\) = 0.5)

Fit the models

extinction_model_sm <-
  brm(data = extinction_data_wrangled,
      family = weibull,
      bf(Gens_to_extinct | cens(Censored_2) ~ 1 + Evolution_treatment + Block + (1|Cross),
         shape ~ 1 + Evolution_treatment + Block),
      # this intercept prior predicts an appropriate starting point for per gen extinction rate on the log scale
      prior = c(prior(gamma(5, 3), class = Intercept, lb = 0),  
                prior(normal(0, 1), class = b),
                prior(exponential(1), class = sd),
                prior(normal(0, 0.5), class = b, dpar = shape)), 
      iter = 8000, warmup = 4000, chains = 4, cores = 4,
      seed = 1, control = list(adapt_delta = .9, max_treedepth = 10),
       file = "Fits/extinction_model")

extinction_model_sm <- add_criterion(extinction_model_sm, criterion = "loo")

# the model does not converge when we try and fit a treatment * block interaction

extinction_model_mm <-
  brm(data = extinction_data_wrangled,
      family = weibull,
      bf(Gens_to_extinct | cens(Censored_2) ~ 1 + Evolution_treatment + Block + (1|mm(Mother_strain, Father_strain)),
      shape ~ 1 + Evolution_treatment + Block),
      # this intercept prior predicts an appropriate starting point for per gen extinction rate on the log scale
      prior = c(prior(gamma(5, 3), class = Intercept, lb = 0), 
                prior(normal(0, 1), class = b),
                prior(exponential(1), class = sd),
                prior(normal(0, 0.5), class = b, dpar = shape)), 
      iter = 8000, warmup = 4000, chains = 4, cores = 4,
      seed = 2, control = list(adapt_delta = 0.98, max_treedepth = 15),
      file = "Fits/extinction_model_mm")

extinction_model_mm <- add_criterion(extinction_model_mm, criterion = "loo")

\(~\)

Model diagnostics

\(~\)

Compare models using LOO, a bayesian information criteria

loo_compare(extinction_model_sm, extinction_model_mm)  %>% 
  kable(digits = 3) %>% 
  kable_styling()

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic
extinction_model_sm	0.000	0.000	-1042.910	20.209	23.694	2.970	2085.820	40.418
extinction_model_mm	-0.668	2.593	-1043.578	20.238	15.072	2.273	2087.155	40.475

The single-membership model performs slightly better. We present results from this model.

Lets see how the model recapitulates the data

pp_check(extinction_model_sm, type = "hist", ndraws = 11, binwidth = 1) +
  theme_minimal() +
  theme(panel.background = element_blank())

\(~\)

Interpreting model output

\(~\)

We can initially ignore the shape parameter and find the mean rate of extinction, averaged across all generations, using the equation

\(\lambda= \frac{1}{e^\mu}\)

First find distributions of \(\mu\) for each of the three treatments.

inverse_rate_weibull <-
  extinction_model_sm %>% 
  as_draws_df() %>% 
  transmute(`Control B1` = b_Intercept, 
         `Female_limited B1` = b_Intercept + b_Evolution_treatmentFemale,
         `Male_limited B1` = b_Intercept + b_Evolution_treatmentMale,
         `Control B2` = b_Intercept + b_Block2, 
         `Female_limited B2` = b_Intercept + b_Evolution_treatmentFemale + b_Block2,
         `Male_limited B2` = b_Intercept + b_Evolution_treatmentMale + b_Block2)

Now lets plug these distributions into the equation to find mean estimates of \(\lambda\).

  inverse_rate_weibull %>% 
  mutate(across(1:6, ~1 / exp(.x))) %>% 
  mutate(across(everything(), ~mean(.x))) %>% 
  distinct()

## # A tibble: 1 × 6
##   `Control B1` `Female_limited B1` `Male_limited B1` Control B…¹ Femal…² Male_…³
##          <dbl>               <dbl>             <dbl>       <dbl>   <dbl>   <dbl>
## 1        0.142               0.130             0.111       0.120   0.110  0.0941
## # … with abbreviated variable names ¹​`Control B2`, ²​`Female_limited B2`,
## #   ³​`Male_limited B2`

This is the rate of extinction per generation, averaged across all generations. However, the weibull family allows a dynamic rate rather the the constant rate that the exponential constrains it to.

Following the vignette("brms_families"), the weibull survival function can be found using the following transformation:

\(\lambda= \frac{e^\mu}{\gamma(1 + \frac{1}{k})}\)

where \(k\) is the shape parameter estimated in brms

Then we can find survival at generation \(t\) using

\(S = e^{-(\frac{t}{\lambda})^k}\)

We plug in a vector with 16000 draws from the posterior distribution for 1) \(\mu\) of each inheritance treatment and 2) \(k\) the shape parameter to find distributions of the proportion of surviving subpopulations at each generation.

# Find lambda for the three treatments. 

lambda_weibull <-
  extinction_model_sm %>% 
  as_draws_df() %>% 
  transmute(`Control B1` = b_Intercept, 
            `Female_limited B1` = b_Intercept + b_Evolution_treatmentFemale,
            `Male_limited B1` = b_Intercept + b_Evolution_treatmentMale,
            `Control B2` = b_Intercept + b_Block2, 
            `Female_limited B2` = b_Intercept + b_Evolution_treatmentFemale + b_Block2,
            `Male_limited B2` = b_Intercept + b_Evolution_treatmentMale + b_Block2,
            shape_control_B1 = exp(b_shape_Intercept),
            shape_female_B1 = exp(b_shape_Intercept + b_shape_Evolution_treatmentFemale),
            shape_male_B1 = exp(b_shape_Intercept + b_shape_Evolution_treatmentMale),
            shape_control_B2 = exp(b_shape_Intercept + b_shape_Block2),
            shape_female_B2 = exp(b_shape_Intercept + b_shape_Evolution_treatmentFemale + b_shape_Block2),
            shape_male_B2 = exp(b_shape_Intercept + b_shape_Evolution_treatmentMale + b_shape_Block2)) %>% 
  mutate(Control_lambda_B1 = exp(`Control B1`) / gamma(1 + 1/shape_control_B1),
         Female_lambda_B1 = exp(`Female_limited B1`) / gamma(1 + 1/shape_female_B1),
         Male_lambda_B1 = exp(`Male_limited B1`) / gamma(1 + 1/shape_male_B1),
         Control_lambda_B2 = exp(`Control B2`) / gamma(1 + 1/shape_control_B2),
         Female_lambda_B2 = exp(`Female_limited B2`) / gamma(1 + 1/shape_female_B2),
         Male_lambda_B2 = exp(`Male_limited B2`) / gamma(1 + 1/shape_male_B2)) %>% 
  select(contains(c("lambda", "shape")))

# Now find survival proportion at generations 1-29

Weibull_survival_curve <- 
  lambda_weibull %>% 
  mutate(#Block 1
    Control_0_B1  = exp(-(0/Control_lambda_B1)^shape_control_B1),  
    Control_1_B1  = exp(-(1/Control_lambda_B1)^shape_control_B1),
    Control_2_B1  = exp(-(2/Control_lambda_B1)^shape_control_B1),
    Control_3_B1  = exp(-(3/Control_lambda_B1)^shape_control_B1),
    Control_4_B1  = exp(-(4/Control_lambda_B1)^shape_control_B1),
    Control_5_B1  = exp(-(5/Control_lambda_B1)^shape_control_B1),
    Control_6_B1  = exp(-(6/Control_lambda_B1)^shape_control_B1),
    Control_7_B1  = exp(-(7/Control_lambda_B1)^shape_control_B1),
    Control_8_B1 = exp(-(8/Control_lambda_B1)^shape_control_B1),
    Control_9_B1 = exp(-(9/Control_lambda_B1)^shape_control_B1),
    Control_10_B1 = exp(-(10/Control_lambda_B1)^shape_control_B1),
    Control_11_B1 = exp(-(11/Control_lambda_B1)^shape_control_B1),
    Control_12_B1 = exp(-(12/Control_lambda_B1)^shape_control_B1),
    Control_13_B1 = exp(-(13/Control_lambda_B1)^shape_control_B1),
    Control_14_B1 = exp(-(14/Control_lambda_B1)^shape_control_B1),
    Control_15_B1 = exp(-(15/Control_lambda_B1)^shape_control_B1),
    Control_16_B1 = exp(-(16/Control_lambda_B1)^shape_control_B1),
    Control_17_B1 = exp(-(17/Control_lambda_B1)^shape_control_B1),
    Control_18_B1 = exp(-(18/Control_lambda_B1)^shape_control_B1),
    Control_19_B1 = exp(-(19/Control_lambda_B1)^shape_control_B1),
    Control_20_B1 = exp(-(20/Control_lambda_B1)^shape_control_B1),
    Control_21_B1 = exp(-(21/Control_lambda_B1)^shape_control_B1),
    Control_22_B1 = exp(-(22/Control_lambda_B1)^shape_control_B1),
    Control_23_B1 = exp(-(23/Control_lambda_B1)^shape_control_B1),
    Control_24_B1 = exp(-(24/Control_lambda_B1)^shape_control_B1),
    Control_25_B1 = exp(-(25/Control_lambda_B1)^shape_control_B1),
    Control_26_B1 = exp(-(26/Control_lambda_B1)^shape_control_B1),
    Control_27_B1 = exp(-(27/Control_lambda_B1)^shape_control_B1),
    Control_28_B1 = exp(-(28/Control_lambda_B1)^shape_control_B1),
    Control_29_B1 = exp(-(29/Control_lambda_B1)^shape_control_B1),
    Female_0_B1  = exp(-(0/Female_lambda_B1)^shape_female_B1),
    Female_1_B1  = exp(-(1/Female_lambda_B1)^shape_female_B1),
    Female_2_B1  = exp(-(2/Female_lambda_B1)^shape_female_B1),
    Female_3_B1  = exp(-(3/Female_lambda_B1)^shape_female_B1),
    Female_4_B1  = exp(-(4/Female_lambda_B1)^shape_female_B1),
    Female_5_B1  = exp(-(5/Female_lambda_B1)^shape_female_B1),
    Female_6_B1  = exp(-(6/Female_lambda_B1)^shape_female_B1),
    Female_7_B1  = exp(-(7/Female_lambda_B1)^shape_female_B1),
    Female_8_B1  = exp(-(8/Female_lambda_B1)^shape_female_B1),
    Female_9_B1  = exp(-(9/Female_lambda_B1)^shape_female_B1),
    Female_10_B1  = exp(-(10/Female_lambda_B1)^shape_female_B1),
    Female_11_B1  = exp(-(11/Female_lambda_B1)^shape_female_B1),
    Female_12_B1  = exp(-(12/Female_lambda_B1)^shape_female_B1),
    Female_13_B1  = exp(-(13/Female_lambda_B1)^shape_female_B1),
    Female_14_B1  = exp(-(14/Female_lambda_B1)^shape_female_B1),
    Female_15_B1  = exp(-(15/Female_lambda_B1)^shape_female_B1),
    Female_16_B1  = exp(-(16/Female_lambda_B1)^shape_female_B1),
    Female_17_B1  = exp(-(17/Female_lambda_B1)^shape_female_B1),
    Female_18_B1  = exp(-(18/Female_lambda_B1)^shape_female_B1),
    Female_19_B1  = exp(-(19/Female_lambda_B1)^shape_female_B1),
    Female_20_B1  = exp(-(20/Female_lambda_B1)^shape_female_B1),
    Female_21_B1 = exp(-(21/Female_lambda_B1)^shape_female_B1),
    Female_22_B1 = exp(-(22/Female_lambda_B1)^shape_female_B1),
    Female_23_B1 = exp(-(23/Female_lambda_B1)^shape_female_B1),
    Female_24_B1 = exp(-(24/Female_lambda_B1)^shape_female_B1),
    Female_25_B1 = exp(-(25/Female_lambda_B1)^shape_female_B1),
    Female_26_B1 = exp(-(26/Female_lambda_B1)^shape_female_B1),
    Female_27_B1 = exp(-(27/Female_lambda_B1)^shape_female_B1),
    Female_28_B1 = exp(-(28/Female_lambda_B1)^shape_female_B1),
    Female_29_B1 = exp(-(29/Female_lambda_B1)^shape_female_B1),
    Male_0_B1     = exp(-(0/Male_lambda_B1)^shape_male_B1),
    Male_1_B1     = exp(-(1/Male_lambda_B1)^shape_male_B1),
    Male_2_B1     = exp(-(2/Male_lambda_B1)^shape_male_B1),
    Male_3_B1     = exp(-(3/Male_lambda_B1)^shape_male_B1),
    Male_4_B1     = exp(-(4/Male_lambda_B1)^shape_male_B1),
    Male_5_B1     = exp(-(5/Male_lambda_B1)^shape_male_B1),
    Male_6_B1     = exp(-(6/Male_lambda_B1)^shape_male_B1),
    Male_7_B1     = exp(-(7/Male_lambda_B1)^shape_male_B1),
    Male_8_B1     = exp(-(8/Male_lambda_B1)^shape_male_B1),
    Male_9_B1     = exp(-(9/Male_lambda_B1)^shape_male_B1),
    Male_10_B1    = exp(-(10/Male_lambda_B1)^shape_male_B1),
    Male_11_B1    = exp(-(11/Male_lambda_B1)^shape_male_B1),
    Male_12_B1    = exp(-(12/Male_lambda_B1)^shape_male_B1),
    Male_13_B1    = exp(-(13/Male_lambda_B1)^shape_male_B1),
    Male_14_B1    = exp(-(14/Male_lambda_B1)^shape_male_B1),
    Male_15_B1    = exp(-(15/Male_lambda_B1)^shape_male_B1),
    Male_16_B1    = exp(-(16/Male_lambda_B1)^shape_male_B1),
    Male_17_B1    = exp(-(17/Male_lambda_B1)^shape_male_B1),
    Male_18_B1    = exp(-(18/Male_lambda_B1)^shape_male_B1),
    Male_19_B1    = exp(-(19/Male_lambda_B1)^shape_male_B1),
    Male_20_B1    = exp(-(20/Male_lambda_B1)^shape_male_B1),
    Male_21_B1 = exp(-(21/Male_lambda_B1)^shape_male_B1),
    Male_22_B1 = exp(-(22/Male_lambda_B1)^shape_male_B1),
    Male_23_B1 = exp(-(23/Male_lambda_B1)^shape_male_B1),
    Male_24_B1 = exp(-(24/Male_lambda_B1)^shape_male_B1),
    Male_25_B1 = exp(-(25/Male_lambda_B1)^shape_male_B1),
    Male_26_B1 = exp(-(26/Male_lambda_B1)^shape_male_B1),
    Male_27_B1 = exp(-(27/Male_lambda_B1)^shape_male_B1),
    Male_28_B1 = exp(-(28/Male_lambda_B1)^shape_male_B1),
    Male_29_B1 = exp(-(29/Male_lambda_B1)^shape_male_B1),
    # Block 2
    Control_0_B2  = exp(-(0/Control_lambda_B2)^shape_control_B2),  
    Control_1_B2  = exp(-(1/Control_lambda_B2)^shape_control_B2),
    Control_2_B2  = exp(-(2/Control_lambda_B2)^shape_control_B2),
    Control_3_B2  = exp(-(3/Control_lambda_B2)^shape_control_B2),
    Control_4_B2  = exp(-(4/Control_lambda_B2)^shape_control_B2),
    Control_5_B2  = exp(-(5/Control_lambda_B2)^shape_control_B2),
    Control_6_B2  = exp(-(6/Control_lambda_B2)^shape_control_B2),
    Control_7_B2  = exp(-(7/Control_lambda_B2)^shape_control_B2),
    Control_8_B2 = exp(-(8/Control_lambda_B2)^shape_control_B2),
    Control_9_B2 = exp(-(9/Control_lambda_B2)^shape_control_B2),
    Control_10_B2 = exp(-(10/Control_lambda_B2)^shape_control_B2),
    Control_11_B2 = exp(-(11/Control_lambda_B2)^shape_control_B2),
    Control_12_B2 = exp(-(12/Control_lambda_B2)^shape_control_B2),
    Control_13_B2 = exp(-(13/Control_lambda_B2)^shape_control_B2),
    Control_14_B2 = exp(-(14/Control_lambda_B2)^shape_control_B2),
    Control_15_B2 = exp(-(15/Control_lambda_B2)^shape_control_B2),
    Control_16_B2 = exp(-(16/Control_lambda_B2)^shape_control_B2),
    Control_17_B2 = exp(-(17/Control_lambda_B2)^shape_control_B2),
    Control_18_B2 = exp(-(18/Control_lambda_B2)^shape_control_B2),
    Control_19_B2 = exp(-(19/Control_lambda_B2)^shape_control_B2),
    Control_20_B2 = exp(-(20/Control_lambda_B2)^shape_control_B2),
    Control_21_B2 = exp(-(21/Control_lambda_B2)^shape_control_B2),
    Control_22_B2 = exp(-(22/Control_lambda_B2)^shape_control_B2),
    Control_23_B2 = exp(-(23/Control_lambda_B2)^shape_control_B2),
    Control_24_B2 = exp(-(24/Control_lambda_B2)^shape_control_B2),
    Control_25_B2 = exp(-(25/Control_lambda_B2)^shape_control_B2),
    Control_26_B2 = exp(-(26/Control_lambda_B2)^shape_control_B2),
    Control_27_B2 = exp(-(27/Control_lambda_B2)^shape_control_B2),
    Control_28_B2 = exp(-(28/Control_lambda_B2)^shape_control_B2),
    Control_29_B2 = exp(-(29/Control_lambda_B2)^shape_control_B2),
    Female_0_B2  = exp(-(0/Female_lambda_B2)^shape_female_B2),
    Female_1_B2  = exp(-(1/Female_lambda_B2)^shape_female_B2),
    Female_2_B2  = exp(-(2/Female_lambda_B2)^shape_female_B2),
    Female_3_B2  = exp(-(3/Female_lambda_B2)^shape_female_B2),
    Female_4_B2  = exp(-(4/Female_lambda_B2)^shape_female_B2),
    Female_5_B2  = exp(-(5/Female_lambda_B2)^shape_female_B2),
    Female_6_B2  = exp(-(6/Female_lambda_B2)^shape_female_B2),
    Female_7_B2  = exp(-(7/Female_lambda_B2)^shape_female_B2),
    Female_8_B2  = exp(-(8/Female_lambda_B2)^shape_female_B2),
    Female_9_B2  = exp(-(9/Female_lambda_B2)^shape_female_B2),
    Female_10_B2  = exp(-(10/Female_lambda_B2)^shape_female_B2),
    Female_11_B2  = exp(-(11/Female_lambda_B2)^shape_female_B2),
    Female_12_B2  = exp(-(12/Female_lambda_B2)^shape_female_B2),
    Female_13_B2  = exp(-(13/Female_lambda_B2)^shape_female_B2),
    Female_14_B2  = exp(-(14/Female_lambda_B2)^shape_female_B2),
    Female_15_B2  = exp(-(15/Female_lambda_B2)^shape_female_B2),
    Female_16_B2  = exp(-(16/Female_lambda_B2)^shape_female_B2),
    Female_17_B2  = exp(-(17/Female_lambda_B2)^shape_female_B2),
    Female_18_B2  = exp(-(18/Female_lambda_B2)^shape_female_B2),
    Female_19_B2  = exp(-(19/Female_lambda_B2)^shape_female_B2),
    Female_20_B2  = exp(-(20/Female_lambda_B2)^shape_female_B2),
    Female_21_B2 = exp(-(21/Female_lambda_B2)^shape_female_B2),
    Female_22_B2 = exp(-(22/Female_lambda_B2)^shape_female_B2),
    Female_23_B2 = exp(-(23/Female_lambda_B2)^shape_female_B2),
    Female_24_B2 = exp(-(24/Female_lambda_B2)^shape_female_B2),
    Female_25_B2 = exp(-(25/Female_lambda_B2)^shape_female_B2),
    Female_26_B2 = exp(-(26/Female_lambda_B2)^shape_female_B2),
    Female_27_B2 = exp(-(27/Female_lambda_B2)^shape_female_B2),
    Female_28_B2 = exp(-(28/Female_lambda_B2)^shape_female_B2),
    Female_29_B2 = exp(-(29/Female_lambda_B2)^shape_female_B2),
    Male_0_B2     = exp(-(0/Male_lambda_B2)^shape_male_B2),
    Male_1_B2     = exp(-(1/Male_lambda_B2)^shape_male_B2),
    Male_2_B2     = exp(-(2/Male_lambda_B2)^shape_male_B2),
    Male_3_B2     = exp(-(3/Male_lambda_B2)^shape_male_B2),
    Male_4_B2     = exp(-(4/Male_lambda_B2)^shape_male_B2),
    Male_5_B2     = exp(-(5/Male_lambda_B2)^shape_male_B2),
    Male_6_B2     = exp(-(6/Male_lambda_B2)^shape_male_B2),
    Male_7_B2     = exp(-(7/Male_lambda_B2)^shape_male_B2),
    Male_8_B2     = exp(-(8/Male_lambda_B2)^shape_male_B2),
    Male_9_B2     = exp(-(9/Male_lambda_B2)^shape_male_B2),
    Male_10_B2    = exp(-(10/Male_lambda_B2)^shape_male_B2),
    Male_11_B2    = exp(-(11/Male_lambda_B2)^shape_male_B2),
    Male_12_B2    = exp(-(12/Male_lambda_B2)^shape_male_B2),
    Male_13_B2    = exp(-(13/Male_lambda_B2)^shape_male_B2),
    Male_14_B2    = exp(-(14/Male_lambda_B2)^shape_male_B2),
    Male_15_B2    = exp(-(15/Male_lambda_B2)^shape_male_B2),
    Male_16_B2    = exp(-(16/Male_lambda_B2)^shape_male_B2),
    Male_17_B2    = exp(-(17/Male_lambda_B2)^shape_male_B2),
    Male_18_B2    = exp(-(18/Male_lambda_B2)^shape_male_B2),
    Male_19_B2    = exp(-(19/Male_lambda_B2)^shape_male_B2),
    Male_20_B2    = exp(-(20/Male_lambda_B2)^shape_male_B2),
    Male_21_B2 = exp(-(21/Male_lambda_B2)^shape_male_B2),
    Male_22_B2 = exp(-(22/Male_lambda_B2)^shape_male_B2),
    Male_23_B2 = exp(-(23/Male_lambda_B2)^shape_male_B2),
    Male_24_B2 = exp(-(24/Male_lambda_B2)^shape_male_B2),
    Male_25_B2 = exp(-(25/Male_lambda_B2)^shape_male_B2),
    Male_26_B2 = exp(-(26/Male_lambda_B2)^shape_male_B2),
    Male_27_B2 = exp(-(27/Male_lambda_B2)^shape_male_B2),
    Male_28_B2 = exp(-(28/Male_lambda_B2)^shape_male_B2),
    Male_29_B2 = exp(-(29/Male_lambda_B2)^shape_male_B2)) %>% 
  select(-c(contains(c("shape", "lambda"))))

Weibull_extinction_estimates_long <-
  Weibull_survival_curve %>% 
  pivot_longer(cols = everything(), names_to = "Evolution treatment", values_to = "Prop_subpopulations_surviving") %>% 
  separate(`Evolution treatment`, into = c("Evolution treatment", "Generation", "Block"), sep = "_") %>% 
    mutate(`Evolution treatment` = case_when(
    `Evolution treatment` == "Female" ~ "Female-limited",
    `Evolution treatment` == "Male" ~ "Male-limited",
    `Evolution treatment` == "Control" ~ "Control"
  ))

\(~\)

Create Figure 1

Make panel a of Figure 1

# values from Falconer and Mackay 1996 pg 90

Inbreeding_coefs <-
  tibble(Generation = 0:20,
         Probability_homozygous = c(0, 0.25, 0.375, 0.5, 0.594, 0.672, 0.734, 0.785, 0.826,
                                    0.859, 0.886, 0.908, 0.926, 0.94, 0.951, 0.961, 0.968, 0.974,
                                    0.979, 0.983, 0.986),
         Probability_pop_fixation = c(0, 0, 0.063, 0.172, 0.293, 0.409, 0.512, 0.601, 0.675, 0.736, 0.785, 
                                      0.826, 0.859, 0.886, 0.908, 0.925, 0.94, 0.951, 0.96, 0.968, 0.975)) %>% 
  mutate(single_locus_heterozygosity = 1 - Probability_homozygous,
         single_locus_not_fixed = 1 - Probability_pop_fixation) 

# for the plot annotations

arrows <- 
  tibble(
    x1 = 3.2,
    x2 = 5.6,
    y1 = 0.28, 
    y2 = 0.49)

# Block 1

weibull_surv_plot_B1 <- 
  Weibull_extinction_estimates_long %>% 
  filter(Block == "B1") %>% 
  group_by(`Evolution treatment`, Generation) %>% 
  tidybayes::median_qi(Prop_subpopulations_surviving, .width = 0.5) %>% 
  mutate(Generation = as.numeric(Generation)) %>% 
  arrange(Generation) %>% 
  # plot!
  ggplot(aes(x = Generation)) +
  geom_line(data = Inbreeding_coefs, aes(y = single_locus_not_fixed), linetype = 3, 
            linewidth = 2, alpha = 1) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper, fill = `Evolution treatment`),
              alpha = 1/2) +
  geom_line(aes(y = Prop_subpopulations_surviving, color = `Evolution treatment`), linetype =5, linewidth = 0.8) +
  
  annotate("text", x = 3.2, y = 0.26, size = 4, label = "Single locus prob. of heterozygosity") +
  geom_curve(
    data = arrows, aes(x = x1, y = y1, xend = x2, yend = y2),
    arrow = arrow(length = unit(0.08, "inch")), size = 0.75,
    color = "gray20", curvature = -0.3) +
  
  scale_fill_manual(values = met.brewer("Hiroshige", 3)) +
  scale_color_manual(values = met.brewer("Hiroshige", 3), guide = "none") +
  scale_x_continuous(breaks = c(0, 5, 10, 15, 20, 25, 30), expand = expansion(mult = c(0.01, 0.01)), limits = c(0, 20)) +
  scale_y_continuous(breaks = c(0, 0.25, .5, 0.75, 1), limits = 0:1, expand = expansion(mult = c(0.01, 0.01))) +
  labs(x = "Generations of inbreeding", y = "Proportion of surviving subpopulations", fill = "Inheritance treatment") +
  theme_bw() + 
  theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
        plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
        legend.background = element_rect(fill='transparent'), #transparent legend bg
        legend.box.background = element_rect(fill='transparent', colour = 'transparent'), #transparent legend panel
        legend.position = c(.95, .95),
        legend.justification = c("right", "top"),
        legend.box.just = "right",
        legend.margin = margin(6, 6, 6, 6),
        text = element_text(size = 14),
        legend.title = element_text(size = 12),
        legend.text = element_text(size = 11))

We can also plot the mean generations until extinction for each evolution treatment. Some data wrangling is first needed.

new_data <- expand_grid(
  Evolution_treatment = extinction_data_wrangled$Evolution_treatment,
  Block = 1:2) %>% 
  distinct(Evolution_treatment, Block) %>%
  mutate(key = paste("V", 1:n(), sep = ""))

weibull_estimates <- 
  fitted(extinction_model_sm, newdata = new_data,
         re_formula = NA, summary = F) %>%  #robust = T) %>% 
  as_tibble() %>% 
  rename(`Female-limited_1` = V1, `Female-limited_2` = V2, `Male-limited_1` = V3, `Male-limited_2` = V4, Control_1 = V5, Control_2 = V6)

mean_weibull_estimates <- 
  weibull_estimates %>% 
  pivot_longer(cols = 1:6, names_to = "Evolution treatment", values_to = "generations") %>% 
  separate(col = `Evolution treatment`, into = c("Evolution treatment", "Block"), sep = "_")


# calculate contrasts

weibull_contrasts_B1 <-
  weibull_estimates %>% 
  mutate(`Female - Control` = `Female-limited_1` - Control_1,
         `Male - Female` = `Male-limited_1` - `Female-limited_1`,
         `Male - Control` = `Male-limited_1` - Control_1) %>% 
  select(-c(`Male-limited_1`, Control_1, `Female-limited_1`,
            `Male-limited_2`, Control_2, `Female-limited_2`)) %>% 
  pivot_longer(cols = 1:3, names_to = "Contrast", values_to = "generation_diff")

weibull_contrasts_B2 <-
  weibull_estimates %>% 
  mutate(`Female - Control` = `Female-limited_2` - Control_2,
         `Male - Female` = `Male-limited_2` - `Female-limited_2`,
         `Male - Control` = `Male-limited_2` - Control_2) %>% 
  select(-c(`Male-limited_1`, Control_1, `Female-limited_1`,
            `Male-limited_2`, Control_2, `Female-limited_2`)) %>% 
  pivot_longer(cols = 1:3, names_to = "Contrast", values_to = "generation_diff")

Create panels b and c of Figure 1

# plot the means

# block 1

ext_p1_B1 <- 
  mean_weibull_estimates %>% 
  filter(Block == "1") %>% 
  ggplot(aes(x = generations, y = `Evolution treatment`)) +
  stat_halfeye(aes(fill = `Evolution treatment`), .width = c(0.66, 0.95), alpha = 1,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.5) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals + # width indicates the uncertainty intervals: here we have 66% and 95% intervals+
  scale_fill_manual(values = met.brewer("Hiroshige", 5)) +
  labs(x = "Generations till extinction", y = "Inheritance treatment") +
  scale_x_continuous(breaks=seq(4, 12, 1), limits = c(4, 12)) +
  theme_bw() + 
  theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
        panel.grid.minor.x = element_blank(), 
        plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
        legend.position = "none", #transparent legend panel
        text = element_text(size=14))

# plot the contrasts

ext_p2_B1 <- 
  weibull_contrasts_B1 %>% 
  ggplot(aes(x = generation_diff, y = fct_relevel(Contrast, "Female - Control", "Male - Female", "Male - Control"))) +
  stat_halfeye(aes(fill = Contrast), .width = c(0.66, 0.95), alpha = 1,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.5) + 
  scale_fill_manual(values = c(carto_pal(7, "Peach")[2], carto_pal(7, "Purp")[1], carto_pal(7, "TealGrn")[1])) +
  geom_vline(xintercept = 0, linetype = 2, colour = "black", size = 1) +
  labs(x = "Diff. in generations till extinction", y = "Treatment contrast") +
  scale_x_continuous(breaks=seq(-4, 6, 2)) +
  theme_bw() + 
  theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
        panel.grid.minor.x = element_blank(),
        plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
        legend.position = "none", #transparent legend panel
        text = element_text(size=14))

We can also estimate the difference between evolutionary treatments in the proportion of subpopulations extant in each generation. Below we calculate each contrast and make panels d, e and f

# Block 1

Difference_gens_B1 <-
  Weibull_survival_curve %>% 
  select(contains("B1")) %>% 
  mutate(G1_diff.c = Male_1_B1 - Control_1_B1,
         G1_diff.f = Male_1_B1 - Female_1_B1,
         G1_diff.fc = Female_1_B1 - Control_1_B1,
         G2_diff.c = Male_2_B1 - Control_2_B1,
         G2_diff.f = Male_2_B1 - Female_2_B1,
         G2_diff.fc = Female_2_B1 - Control_2_B1,
         G3_diff.c = Male_3_B1 - Control_3_B1,
         G3_diff.f = Male_3_B1 - Female_3_B1,
         G3_diff.fc = Female_3_B1 - Control_3_B1,
         G4_diff.c = Male_4_B1 - Control_4_B1,
         G4_diff.f = Male_4_B1 - Female_4_B1,
         G4_diff.fc = Female_4_B1 - Control_4_B1,
         G5_diff.c = Male_5_B1 - Control_5_B1,
         G5_diff.f = Male_5_B1 - Female_5_B1,
         G5_diff.fc = Female_5_B1 - Control_5_B1,
         G6_diff.c = Male_6_B1 - Control_6_B1,
         G6_diff.f = Male_6_B1 - Female_6_B1,
         G6_diff.fc = Female_6_B1 - Control_6_B1,
         G7_diff.c = Male_7_B1 - Control_7_B1,
         G7_diff.f = Male_7_B1 - Female_7_B1,
         G7_diff.fc = Female_7_B1 - Control_7_B1,
         G8_diff.c = Male_8_B1 - Control_8_B1,
         G8_diff.f = Male_8_B1 - Female_8_B1,
         G8_diff.fc = Female_8_B1 - Control_8_B1,
         G9_diff.c = Male_9_B1 - Control_9_B1,
         G9_diff.f = Male_9_B1 - Female_9_B1,
         G9_diff.fc = Female_9_B1 - Control_9_B1,
         G10_diff.c = Male_10_B1 - Control_10_B1,
         G10_diff.f = Male_10_B1 - Female_10_B1,
         G10_diff.fc = Female_10_B1 - Control_10_B1,
         G11_diff.c = Male_11_B1 - Control_11_B1,
         G11_diff.f = Male_11_B1 - Female_11_B1,
         G11_diff.fc = Female_11_B1 - Control_11_B1,
         G12_diff.c = Male_12_B1 - Control_12_B1,
         G12_diff.f = Male_12_B1 - Female_12_B1,
         G12_diff.fc = Female_12_B1 - Control_12_B1,
         G13_diff.c = Male_13_B1 - Control_13_B1,
         G13_diff.f = Male_13_B1 - Female_13_B1,
         G13_diff.fc = Female_13_B1 - Control_13_B1,
         G14_diff.c = Male_14_B1 - Control_14_B1,
         G14_diff.f = Male_14_B1 - Female_14_B1,
         G14_diff.fc = Female_14_B1 - Control_14_B1,
         G15_diff.c = Male_15_B1 - Control_15_B1,
         G15_diff.f = Male_15_B1 - Female_15_B1,
         G15_diff.fc = Female_15_B1 - Control_15_B1,
         G16_diff.c = Male_16_B1 - Control_16_B1,
         G16_diff.f = Male_16_B1 - Female_16_B1,
         G16_diff.fc = Female_16_B1 - Control_16_B1,
         G17_diff.c = Male_17_B1 - Control_17_B1,
         G17_diff.f = Male_17_B1 - Female_17_B1,
         G17_diff.fc = Female_17_B1 - Control_17_B1,
         G18_diff.c = Male_18_B1 - Control_18_B1,
         G18_diff.f = Male_18_B1 - Female_18_B1,
         G18_diff.fc = Female_18_B1 - Control_18_B1,
         G19_diff.c = Male_19_B1 - Control_19_B1,
         G19_diff.f = Male_19_B1 - Female_19_B1,
         G19_diff.fc = Female_19_B1 - Control_19_B1,
         G20_diff.c = Male_20_B1 - Control_20_B1,
         G20_diff.f = Male_20_B1 - Female_20_B1,
         G20_diff.fc = Female_20_B1 - Control_20_B1) %>% 
        #G21_diff.c = Male_21_B1 - Control_21_B1,
        #G21_diff.f = Male_21_B1 - Female_21_B1,
        #G21_diff.fc = Female_21_B1 - Control_21_B1,
        #G22_diff.c = Male_22_B1 - Control_22_B1,
        #G22_diff.f = Male_22_B1 - Female_22_B1,
        #G22_diff.fc = Female_22_B1 - Control_22_B1,
        #G23_diff.c = Male_23_B1 - Control_23_B1,
        #G23_diff.f = Male_23_B1 - Female_23_B1,
        #G23_diff.fc = Female_23_B1 - Control_23_B1,
        #G24_diff.c = Male_24_B1 - Control_24_B1,
        #G24_diff.f = Male_24_B1 - Female_24_B1,
        #G24_diff.fc = Female_24_B1 - Control_24_B1,
        #G25_diff.c = Male_25_B1 - Control_25_B1,
        #G25_diff.f = Male_25_B1 - Female_25_B1,
        #G25_diff.fc = Female_25_B1 - Control_25_B1,
        #G26_diff.c = Male_26_B1 - Control_26_B1,
        #G26_diff.f = Male_26_B1 - Female_26_B1,
        #G26_diff.fc = Female_26_B1 - Control_26_B1,
        #G27_diff.c = Male_27_B1 - Control_27_B1,
        #G27_diff.f = Male_27_B1 - Female_27_B1,
        #G27_diff.fc = Female_27_B1 - Control_27_B1,
        #G28_diff.c = Male_28_B1 - Control_28_B1,
        #G28_diff.f = Male_28_B1 - Female_28_B1,
        #G28_diff.fc = Female_28_B1 - Control_28_B1,
        #G29_diff.c = Male_2_B1 - Control_29_B1,
        #G29_diff.f = Male_29_B1 - Female_29_B1,
        #G29_diff.fc = Female_29_B1 - Control_29_B1) %>% 
  select(starts_with("G")) %>% 
  pivot_longer(cols = everything(), names_to = "Difference contrast", values_to = "survival_diff") %>% 
  separate(`Difference contrast`, into = c("Generation", "Difference contrast"), sep = "_") %>% 
  mutate(Generation = as.numeric(str_remove(Generation, "G")))


# Make the three plots

Leaf_plot_mc_B1 <-
  Difference_gens_B1 %>%
  filter(`Difference contrast` == "diff.c") %>% 
  ggplot(aes(survival_diff, as.factor(Generation))) +
  stat_interval(.width = c(0.05, 0.66, 0.95), 
                height = 1, show.legend = F) +
  rcartocolor::scale_color_carto_d(palette = "Purp") +
  coord_cartesian(xlim = c(-0.1, 0.3)) +
  scale_y_discrete(limits=rev) +
  geom_vline(linetype = 2, xintercept = 0, size = 1) +
  labs(x = "Male - Control\nProp. surv subpopulations",
       y = "Generation of inbreeding") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size=14),
        axis.text.y = element_text(size = 8)) 

Leaf_plot_mf_B1 <-
  Difference_gens_B1 %>%
  filter(`Difference contrast` == "diff.f") %>% 
  ggplot(aes(survival_diff, as.factor(Generation))) +
  stat_interval(.width = c(0.05, 0.66, 0.95), 
                height = 1, show.legend = F) +
  rcartocolor::scale_color_carto_d(palette = "TealGrn") +
  coord_cartesian(xlim = c(-0.1, 0.3)) +
  scale_y_discrete(limits=rev) +
  geom_vline(linetype = 2, xintercept = 0, size = 1) +
  labs(x = "Male - Female\nProp. surv subpopulations",
       y = NULL) +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size=14),
        axis.text.y = element_text(size = 8))

Leaf_plot_fc_B1 <-
  Difference_gens_B1 %>%
  filter(`Difference contrast` == "diff.fc") %>% 
  ggplot(aes(survival_diff, as.factor(Generation))) +
  stat_interval(.width = c(0.05, 0.66, 0.95), 
                height = 1, show.legend = F) +
  rcartocolor::scale_color_carto_d(palette = "Peach") +
  coord_cartesian(xlim = c(-0.1, 0.3)) +
  scale_y_discrete(limits=rev) +
  geom_vline(linetype = 2, xintercept = 0, size = 1) +
  labs(x = "Female - Control\nProp. surv subpopulations",
       y = NULL) +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size=14),
        axis.text.y = element_text(size = 8))

Combine the 6 panels

(weibull_surv_plot_B1 / (ext_p1_B1 + ext_p2_B1)) / (Leaf_plot_mc_B1 + Leaf_plot_mf_B1 + Leaf_plot_fc_B1) +
  plot_annotation(tag_levels = 'a') + 
  plot_layout(heights = c(1.5, 1, 1))

Figure 1. Panel a shows extinction trajectories under increasing stress from inbreeding. Bands are 50% credible intervals to aid visual inspection. The black dotted line is the probability that a randomly selected autosomal locus remains heterozygous. b shows the posterior distribution for the number of generations that sub-populations descended from each inheritance treatment survived inbreeding exposure. The white point shows the median with associated 66 and 95% credible intervals c shows the distribution of the difference between each treatment in generations survived. d-f show the distribution of the difference between each treatment in generations survived, estimated separately for each generation of inbreeding. The innermost colour band approximates the median, while the middle and outermost bars show 66 and 95% credible intervals.

\(~\)

Productivity analysis

\(~\)

We have several hypotheses that we test with the productivity data. First, productivity should be negatively correlated with the number of generations inbreeding occurs, due to inbreeding depression. As mentioned above, this effect should be strongest in the first ~10 generations of the experiment, during which full-sibling inbreeding should have the largest effects on genome wide heterozygosity. After this point, the inbreeding coefficient will be 0.89, while over the next 10 generations it should only increase by a further 0.1 to 0.99. Hence, we expect productivity to stabilise between the 10th and 20th generations, assuming the subpopulation is still extant. Second, subpopulations carrying autosomes that have responded to selection on males should exhibit a smaller drop in productivity compared with subpopulations carrying the female-adapted or control autosomes, because they have a history of stronger selection that should, in theory, have purged the genome of recessive deleterious alleles.

Load in the data

data <-
  read_csv("Data/Productivity_data.csv") %>% 
  mutate(Week = case_when(
    Block == 1 ~ Generation,
    Block == 2 ~ Generation + 9))

Productivity_data <- 
  left_join(
    data,
    
    data %>% 
      distinct(Cross, Replicate) %>% 
      rowid_to_column("subpopulation")
  ) %>% 
  select(subpopulation, everything()) %>% 
  mutate(across(1:7, as.factor),
         Week = as.factor(Week),
         Collection_window_offspring = Female_offspring + Male_offspring,
         Pre_window_offspring = Pre_window_female_offspring + Pre_window_male_offspring,
         Total_female_offspring = Female_offspring + Pre_window_female_offspring,
         Total_male_offspring = Male_offspring + Pre_window_male_offspring,
         Total_offspring = Total_female_offspring + Total_male_offspring) %>% 
  # In one generation of the experiment (b1 = G24 & B2 = G15) sibling pairs were setup a day early and removed from their vials at the regular time, meaning that they had an extra day to produce offspring. To correct for this we multiply offspring counts by 0.75.
  mutate(Collection_window_offspring = if_else(Count_conditions == "Extra day", 
                                               round(Collection_window_offspring * 0.75),  Collection_window_offspring),
         # note that because everything is moved a day early, pre-window offspring counts will still be inflated even after this correction
         Pre_window_offspring = if_else(Count_conditions == "Extra day", 
                                        round(Pre_window_offspring * 0.75),  Pre_window_offspring))
  
# Combine with mortality and extinction data to filter out vials that went extinct because of female mortality

Productivity_mortality <-
  left_join(
    Productivity_data %>% 
      filter(Generation < 21), #%>% only include data collected on the first 20 generations of inbreeding (block 2's endpoint)
    mortality_data %>% select(subpopulation, Generation, Female_mortality)
  )

# We also include a column that specifies if a subpopulation was extinct. This means we can easily remove 0 values from the data if required. We can use the `extinction_data_wrangled$Gens_to_extinct` column to help us here. 

Productivity_data_clean <-
  left_join(
    Productivity_mortality,  
    extinction_data_wrangled %>% 
      select(subpopulation, Gens_to_extinct, Censored_mortality)
  ) %>% 
  mutate(Vial_setup = if_else(Generation > Gens_to_extinct, "No", "Yes")) %>% 
  filter(Collection_window_offspring != "NA")

# remove mortality events and their consequences in later generations

Productivity_data_clean_2 <-
  Productivity_data_clean %>%
  mutate(Collection_window_offspring = case_when(
    Censored_mortality == 1 & Vial_setup == "No" ~ NA_real_,
    TRUE ~ Collection_window_offspring)) %>% 
  filter(!is.na(Collection_window_offspring),
         Female_mortality != 1 | is.na(Female_mortality))


# Create a function to build HTML searchable tables

my_data_table <- function(df){
  datatable(
    df, rownames=FALSE,
    autoHideNavigation = TRUE,
    extensions = c("Scroller",  "Buttons"),
    options = list(
      dom = 'Bfrtip',
      deferRender=TRUE,
      scrollX=TRUE, scrollY=400,
      scrollCollapse=TRUE,
      buttons =
        list('pageLength', 'colvis', 'csv', list(
          extend = 'pdf',
          pageSize = 'A4',
          orientation = 'landscape',
          filename = 'Productivity_data_clean_2')),
      pageLength = 8500
    )
  )
}

my_data_table(Productivity_data_clean %>% 
                select(Mother_strain, Father_strain, Cross, subpopulation, Block, Week, Generation, Treatment, Vial_setup, Count_conditions, Female_offspring, Male_offspring, Collection_window_offspring, Pre_window_female_offspring, Pre_window_male_offspring, Pre_window_offspring, Total_female_offspring, Total_male_offspring, Total_offspring))

Column explanations

Mother strain: we measured extinction and productivity of families derived from 36 different crosses. Crosses were only conducted between strains that experienced the same Evolution_treatment. This column shows the strain that the founding female of the subpopulation was derived from.

Father strain: we measured extinction and productivity of families derived from 36 different crosses. Crosses were only conducted between strains that experienced the same evolution treatment. This column shows the strain that the founding male of the subpopulation was derived from.

Cross: an identifying ID for each of the 36 combinations of Mother strain and Father strain.

subpopulation: lines of flies founded by single inseminated females in generation zero. We generated the next generation of each subpopulation by crossing a single full-sibling dyad.

Block: the experiment was run in 2 distinct blocks, using flies separated by 9 generations. Block 1 ran for 29 generations, whereas Block 2 ran for 20 generations.

Week: the week that the generation of each block was conducted in. Note that blocks overlap between weeks 9-29.

Generation: the generation of the extinction assay, ranging from 1 to 20.

Treatment: the strains had been exposed to one of three evolutionary conditions for 20-29 generations: a female-limited response to selection, a male-limited response and a control condition where an evolutionary response occurred in both sexes.

Vial setup: specifies whether a vial was setup for a given subpopulation in a given generation. Vials were not setup when the subpopulation was extinct, or when it could no longer be included in the experiment due to escape/ handling error.

Count_conditions: Normal indicates that the breeding window lasted the standard three days. Extra day indicates the breeding pair were left in the vial for an extra day. This occurred during generation 15 of the second experimental block.

Female_offspring: the number of adult females produced by the breeding pair that eclosed during the four-day collection window

Male_offspring: as above, but for male offspring.

Collection_window_offspring: the sum of Female_offspring and Male_offspring.

Pre_window_female_offspring: females progeny that eclosed prior to the collection window opening.

Pre_window_male_offspring: as above, but for male offspring.

Pre_window_offspring: the sum of Pre_window_female_offspring and Pre_window_male_offspring.

Total_female_offspring: the sum of Female_offspring and Pre_window_female_offspring.

Total_male_offspring: as above, but for male offspring.

Total_offspring: the total number of offspring produced, summed across the sexes and the pre-window and within-window collections.

\(~\)

Do offspring emerge prior to day 10 throughout the experiment?

  Productivity_data_clean %>% 
  ggplot(aes(x = Generation, y = Pre_window_offspring/Total_offspring)) +
  geom_jitter(aes(size = Total_offspring), colour = met.brewer("OKeeffe2")[2], fill = met.brewer("OKeeffe2")[3], width = 0.2, 
              shape = 21, stroke =1, alpha = 0.25) +
  scale_size(range = c(0, 5), name = "Total offspring") +
  scale_y_continuous(limits = c(0,1), expand = expansion(mult = c(0.01, 0.01))) +
  theme_tidybayes() +
  coord_cartesian(xlim = c(1, 20)) +
  labs(y = "Prop. offspring eclosing prior to day 10", 
       x = "Generation of inbreeding",
       ) +
  theme(text = element_text(size = 14),
        legend.position = "bottom")

Figure S3. Raw proportion of productivity that occurred prior to the eclosion window, for the first 20 generations of the experiment.

\(~\)

Testing for our casual effect

\(~\)

Modelling approach

As mentioned above, we have strong expectations for how full-sibling inbreeding should affect productivity. Decreasing genome-wide heterozygosity and the expression of recessive deleterious alleles that comes with this should follow a pattern of exponential decay. This non-linear process can be modelled using the brms non-linear syntax.

Here we model the productivity at generation \(t\) via the function \(a - e^{\lambda t}\) where \(a\) and \(\lambda\) are parameters that control the initial productivity and the rate of its decay as generations of inbreeding progress. We force \(a\) to be positive using the exp() function, which expresses it on the exponential scale. Large \(a\) indicates higher starting productivity at generation 0 (we started observation at generation 1, the productivity for which is given by \(a - e^{\lambda}\)). Positive values of \(\lambda\) indicate that productivity declines with inbreeding, which is our strong expectation and clearly the case following inspection of the raw data (see Figure SX).

We model inbreeding depression while including zero values for all rows where a subpopulation was extinct. This helps combat underestimation of the decay parameter, caused by the non-random extinction of subpopulations, where high fitness subpopulations become over-represented as the generations of inbreeding progress (i.e. selection).

Fixed and random effects

In the brms non-linear syntax, different formulas can be fit for each parameter. We fit the same fixed effect for \(a\) and \(\lambda\): Evolution treatment. We fit subpopulation as a random effect to account for repeated measurements taken on each subpopulation for \(a\), but are unable to do this for \(\lambda\), as individual subpopulation productivity curves are highly heteroskedastic, and are skewed by extinction events. To minimise pseudoreplication for \(\lambda\), we fit Cross as a random effect, as we do for the extinction data. We also fit Week as a random effect for both parameters; this controls for environmental variation between generations.

Priors

We fit moderately informative priors, with the aim to regularise our posterior estimates and improve model fitting by ruling out non-sensical values. Note that priors are particularly important when fitting non-linear models.

Note that \(\alpha\) related priors are expressed on the exponential scale.

\(\alpha\) (initial productivity) ~ Normal(\(\mu\) = 4.2, \(\sigma\) = 0.25)

\(\lambda\) (the rate of decay) ~ Gamma(shape = 0.8, rate = 4)

\(\sigma\) for \(\alpha\) ~ Gamma(shape = 2, rate = 1.5)

\(\sigma\) for \(\lambda\) ~ Gamma(shape = 0.25, rate = 2)

\(~\)

Simulate the curve, with the priors we will specify

Prior predictive simulation for the productivity decline model. The plot shows the results from 100 prior draws.

n <- 1e4

# define the x-axis breaks
at <- 1:20

# how many prior draws would you like?
n_draws <- 200

# simulate
set.seed(16)

prior <-
  tibble(index = 1:n,
         a   = rnorm(n, mean = 4.2, sd = 0.25),
         lambda = rgamma(n, 0.8, 4)) %>% 
  slice_sample(n = n_draws) %>% 
  expand(nesting(index, a, lambda),
         Generation = seq(from = 0, to = 20, length.out = 1e2)) %>% 
  mutate(productivity = exp(a - lambda * Generation)) 


prior %>% 
  ggplot(aes(x = Generation, y = productivity, group = index)) +
  geom_line(size = 1/2, alpha = 1/4, colour = met.brewer("Hiroshige")[2]) +
  theme_tidybayes() +
  scale_y_continuous(limits = c(0,200), expand = expansion(mult = c(0, .1))) +
  scale_x_continuous(limits = c(0,20), expand = c(0, 0)) +
  labs(y = "subpopulation productivity", x = "Generation of inbreeding") +
  theme(text = element_text(size = 14))

\(~\)

Fit the non-linear model

nonlinear_productivity <- 
  brm(bf(Collection_window_offspring ~ exp(a - (lambda * Generation)),
         a ~ 1 + Treatment + (1|Week) + (1|subpopulation), 
         lambda ~ 1 + Treatment + (1|Week) + (1|Cross), 
         nl = TRUE),
      data = Productivity_data_clean,
      prior = c(prior(normal(4.2, 0.25), class = b, coef = Intercept, nlpar = "a"),
                prior(normal(0, 0.25), class = b, coef = TreatmentFemale, nlpar = "a"),
                prior(normal(0, 0.25), class = b, coef = TreatmentMale, nlpar = "a"),
                prior(gamma(0.8, 4), class = b, coef = Intercept, nlpar = "lambda"),
                prior(normal(0, 0.2), class = b, coef = TreatmentFemale, nlpar = "lambda"),
                prior(normal(0, 0.2), class = b, coef = TreatmentMale, nlpar = "lambda"),
                prior(gamma(2, 1.5), class = sd, nlpar = "a"),
                prior(gamma(0.25, 2), class = sd, nlpar = "lambda")),
      control = list(adapt_delta = 0.9, max_treedepth = 12), seed = 5,
      iter = 20000, warmup = 5000, chains = 4, cores = 4,
      file = "Fits/non_linear_productivity_model_no_fam")

\(~\)

Build figure 2

\(~\)

new_data <- 
  Productivity_data_clean_2 %>% 
  select(Generation, Treatment) %>%
  distinct() %>%
  mutate(key = paste("V", 1:n(), sep = ""))

preds <- 
  as.data.frame(fitted(nonlinear_productivity, newdata = new_data, re_formula = NA, summary=F)) %>% 
  mutate(draw = 1:n()) %>% 
  gather(key, Collection_window_offspring, -draw) %>% 
  as_tibble() %>% 
  left_join(new_data, by = "key") %>% 
  select(-key) %>% 
    mutate(Treatment = case_when(
    Treatment == "Female" ~ "Female-limited",
    Treatment == "Male" ~ "Male-limited",
    Treatment == "Control" ~ "Control"))

preds_summary <-
  preds %>% 
  group_by(Treatment, Generation) %>%
  median_qi(Collection_window_offspring, .width = 0.5)

nl_p1 <-
  Productivity_data_clean_2 %>% 
  ggplot(aes(x = Generation, y = Collection_window_offspring)) +
  geom_jitter(colour = met.brewer("OKeeffe2")[2], width = 0.25, 
              shape = 1.5, stroke =2, size = 1.2, alpha = 0.2) +
 # geom_smooth(size = 1.5, colour = met.brewer("OKeeffe2")[7], alpha = 0.4) +
  geom_ribbon(data = preds_summary, 
             aes(ymin = .lower, ymax = .upper, fill = Treatment), alpha = 1/2) +
  geom_line(aes(group = Treatment, colour = Treatment), data = preds_summary,
            size = 1, alpha = 1, linetype = 5) +
  scale_linetype_manual(values = c(1, 2, 3)) +
  scale_colour_manual(values = met.brewer("Hiroshige", 3), guide = "none") +
  scale_fill_manual(values = met.brewer("Hiroshige", 3)) +
  scale_y_continuous(expand = expansion(mult = c(0.01, 0))) +
  theme_tidybayes() +
  coord_cartesian(xlim = c(1, 20), ylim = c(0, 150)) +
  labs(y = "Subpopulation productivity", x = "Generation of inbreeding", fill = "Inheritance treatment") +
  theme(text = element_text(size = 14),
        legend.title = element_text(size = 12),
        legend.text = element_text(size = 11),
        legend.position = c(.99, .99),
        legend.justification = c("right", "top"),
        legend.box.just = "right",
        legend.margin = margin(6, 6, 6, 6))

\(~\)

Estimating Generation one productivity and \(\lambda\)

Our model estimates \(a\) the initial productivity parameter. This is the productivity in generation 0, which in practice is non-sensical. We instead present productivity in generation one which = \(a - e^{\lambda*1}\).

Gen_1_productivity <-
preds %>% 
  filter(Generation == 1)
  

Initial_productivity_plot <-
  Gen_1_productivity %>% 
  
  ggplot(aes(x = Collection_window_offspring, y = Treatment)) +
  stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 1,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.5) + 
  scale_fill_manual(values = met.brewer("Hiroshige", 5)) +
  labs(x= "Productivity in generation 1", y = "Inheritance treatment") +
  theme_bw() + 
  theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
        plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
        panel.grid.minor.x = element_blank(),
        legend.position = "none", #transparent legend panel
        text = element_text(size=14))



decay_estimates <-  
  as_draws_df(nonlinear_productivity, variable = "^b_lambda", regex = TRUE) %>% 
  mutate(`Female-limited` = b_lambda_Intercept + b_lambda_TreatmentFemale,
         `Male-limited` = b_lambda_Intercept + b_lambda_TreatmentMale,
         Control = b_lambda_Intercept) %>%
 select(`Female-limited`, `Male-limited`, Control) %>% 
  mutate(posterior_sample = 1:n()) %>% 
  gather(Evolution_treatment, lambda, -posterior_sample)
  

decay_plot <-
  decay_estimates %>% 
  ggplot(aes(x = lambda, y = Evolution_treatment)) +
    stat_halfeye(aes(fill = Evolution_treatment), .width = c(0.66, 0.95), alpha = 1,
               point_interval = "mean_qi", point_fill = "white", 
               shape = 21, point_size = 4, stroke = 1.5, scale = 0.8) +
  scale_fill_manual(values = met.brewer("Hiroshige", 5)) +
  labs(x=expression(paste("Decay rate parameter ", lambda,)), y = NULL) +
  theme_bw() + 
  theme(panel.background = element_rect(fill='transparent'), #transparent panel bg
        plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
        panel.grid.minor.x = element_blank(),
        legend.position = "none", #transparent legend panel
        text = element_text(size=14))

Create plots d and e

Difference_gens_productivity <- 
  preds %>% 
  pivot_wider(values_from = "Collection_window_offspring", names_from = "Treatment") %>% 
  mutate(diff.mc = `Male-limited` - Control,
         diff.mf = `Male-limited` - `Female-limited`,
         diff.fc = `Female-limited` - Control) %>% 
  select(Generation, contains("diff")) %>% 
  pivot_longer(cols = contains("diff"), values_to = "productivity_diff", names_to = "Difference contrast") 

# Make the three plots

Leaf_plot_mc_productivity <-
  Difference_gens_productivity %>%
  filter(`Difference contrast` == "diff.mc") %>% 
  ggplot(aes(productivity_diff, as.factor(Generation))) +
  stat_interval(.width = c(0.05, 0.66, 0.95), 
                height = 1, show.legend = F) +
  rcartocolor::scale_color_carto_d(palette = "Purp") +
  coord_cartesian(xlim = c(-5, 25)) +
  scale_y_discrete(limits=rev) +
  geom_vline(linetype = 2, xintercept = 0, size = 1) +
  labs(x = "Male - Control\nProductivity diff.",
       y = "Generation of inbreeding") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size=14),
        axis.text.y = element_text(size = 8))


Leaf_plot_mf_productivity <-
  Difference_gens_productivity %>%
  filter(`Difference contrast` == "diff.mf") %>% 
  ggplot(aes(productivity_diff, as.factor(Generation))) +
  stat_interval(.width = c(0.05, 0.66, 0.95), 
                height = 1, show.legend = F) +
  rcartocolor::scale_color_carto_d(palette = "TealGrn") +
  coord_cartesian(xlim = c(-5, 25)) +
  scale_y_discrete(limits=rev) +
  geom_vline(linetype = 2, xintercept = 0, size = 1) +
  labs(x = "Male - Female\nProductivity diff.",
       y = NULL) +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size=14),
        axis.text.y = element_text(size = 8))


Leaf_plot_fc_productivity <-
  Difference_gens_productivity %>%
  filter(`Difference contrast` == "diff.fc") %>% 
  ggplot(aes(productivity_diff, as.factor(Generation))) +
  stat_interval(.width = c(0.05, 0.66, 0.95), 
                height = 1, show.legend = F) +
  rcartocolor::scale_color_carto_d(palette = "Peach") +
  coord_cartesian(xlim = c(-5, 25)) +
  scale_y_discrete(limits=rev) +
  geom_vline(linetype = 2, xintercept = 0, size = 1) +
  labs(x = "Female - Control\nProductivity diff.",
       y = NULL) +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size=14),
        axis.text.y = element_text(size = 8))

Combine the 6 panels

(nl_p1 / (Initial_productivity_plot + decay_plot)) / (Leaf_plot_mc_productivity + Leaf_plot_mf_productivity + Leaf_plot_fc_productivity) +
  plot_annotation(tag_levels = 'a') + 
  plot_layout(heights = c(1.5, 1, 1))

Figure 2. Panel a shows productivity decline under increasing stress from inbreeding, split by inheritance treatment. Bands are 50% credible intervals to aid visual inspection and circular points show raw productivity counts for each sub-population in each generation. b shows the posterior distribution of productivity in the first generation of the experiment; white points are the median with associated 66 and 95% credible intervals. c shows the distribution of the productivity decay rate for each treatment. d-f show the difference contrasts for productivity in each generation. The innermost colour band approximates the median, while the middle and outermost bars show 66 and 95% credible intervals.